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Abstract 

This paper treats elliptic problems with corner singularities. Finite 
element approximations based on variational principles of the least squares 
type tend to display poor convergence properties in such contexts. Moreover, 
mesh refinement or the use of special singular elements do not appreciably 
improve matters. Here we show that if the least squares formulation is done 
in appropriately weighted space, then optimal convergence results in 
unweighted spaces like L 2 . 
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!• Introduction 


Least squares methods have proven to be useful for indefinite elliptic 
systems. The Helmholtz equation is perhaps the most important example ( [ 1 ] — 
[2]). The main advantage over standard Galerkin formulations is that least 
squares always yield Hermitian definite algebraic systems. Thus iterative 
methods like SOR can be used [3] . 

The main disadvantage of least squares is the extreme regularity 
requirement it has for optimal rates of convergence. For example, the 
standard least squares approximation to 

(1*1) A<}> + q(J> = f in ft 

(1*2) = 0 on 9ft 

is to require that 

0*3) J { |grad <{> - u | 2 + | div u + q<j> - f | 2 } 

J2 

be minimized as <|> and u vary over appropriate finite dimensional spaces. 
It has been shown that such an approach will give optimal convergence 

(i.e., second-order if linear elements are used, etc.) only in special 
circumstances [1]. The most restrictive condition being the existence of a 
number 0 < c < 09 such that for any f in the Sobolev space H^(£3) the 

solution <)> of ( 1 . 1)— ( 1.2) lies in H 3 (fl) and 

<!•*) H<f>H < c Ilf II 

H J (n) H^G) 
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A regularity result of this type is valid only for smooth regions ft, and 
in particular corners on 3ft are excluded. Numerical experiments suggest 
that this condition may in fact be necessary. For example, a series of 
calculations dealing with rectangular polygons in IF? having re-entrant 
corners showed that not only was the Ir convergence suboptimal on uniform 
grids but it also remained suboptimal even in the presence of mesh refinement 

[4]. 

In this paper we consider an alternative least squares approximation in 
weighted Sobolev spaces. These are spaces where analogs of (1*4) are valid. 
The key feature of our analysis is that the error estimates are in unweighted 
norms like L^. Selected numerical experiments with this type of formulation 
and with appropriate mesh refinement are also reported here. 

The results obtained here generalize those in [5]. As in the latter the 
Hardy-Littlewood inequality plays a key role; however, in this paper the 
analysis takes a different direction in the sense that the discrete 
decomposition property introduced in [6] is also used extensively. In 
addition, the special mesh refinement introduced by Babuska, Kellogg, 
and Pitkaranta [7] is also exploited explicitly. 

For simplicity we shall consider planar regions ft having only one 
corner as is shown in Figure 1-1. Since existing proofs of the Hardy- 
Littlewood inequality [8] uses both exterior and interior cone conditions our 
results are restricted to interior angles 0^ satisfying 


0 < e Q < 2 tt. 
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Figure 1—1: The planar region ft. 


2. Variational Formulation 

Let r denote the distance to the corner at P Q (see Figure 1-1 
9ft. For a nonnegative integer k and a nonnegative number 0 let 


( 2 . 1 ) 


l\pll 


k,0 


k 

- I I 

j=0 ft 


: 2 (S+ j" k) V „,| 2 , 


and let fi k,e (ft) denote the closure of [C (ft)] 2 in this norm. We 

o. 

(J .1)— (1.2) with a least squares formulation in the space H X (ft) 
(where H^(ft) is space of functions in H*(ft) with zero trace on 
particular, let 

< 2 ‘ 2 ) S h CH 1 (ft), \f h CW 1,e (ft) 


) on 


approach 
9ft). In 


be finite dimensional subspaces parameterized by h > 0. We seek 
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( 2.3) V £ s h , u h £ 5 h 

which minimize 

(2.4) / { I grad - vj ^ + r ^ | div v_ + qij; - f | ^ } 

a 

+ 26 
over (^,v_) € x V^. That is, we have (1.3) with a weight r P on the most 

highly differential terms. Our goal is to find appropriate g for which 

9 

<f> ~ converge in unweighted spaces like L at the optimal rates. 

Our analysis can be also used for the case where one has weights on both 
terms in (2.4). It can be shown that the weight on the first term does not 
help, and to minimize technical details we anticipate this result and start 
with (2.4). An intuitive justification for (2.4) can be obtained from the 
nature of the corner singularity. Indeed, if 0^ denotes the interior angle, 
then the solution <j> to (1.1) — (1.2) will have a singularity of the form 

Tr/ 0 ir/e n -l 

(2.5) <J> ~ 0(r ), 11 = grad $ ~ 0(r ) 


and here we are concerned with the case or re-entrant corners where tt < 0^. 
Thus if one were dealing with functions (\(j,v) in S,^ having this type of 
behavior (such as special singular elements), then 


( 2 . 6 ) 


11 / 9 n 

0(r °), 


TT/0 — i 

v = 0(r 0 ) 


and so 


7T/0Q-2 


). 


(2.7) 


div v = 0(r 
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Observe that grad \|> and _v are square integrable; hence the first term in 
(2.4) is well defined. However, div v_ is not square integrable, and the 

second term in (2.4) is finite only if 8 exceeds 1 - ir/0Q. Our analysis 
(Section 4) indicates that if 8 is slightly larger than this; i.e., 

(2.8) 1 - tt / 2 0 Q < 0 < 1; 

then optional rates of convergence will result in L2 under suitable 

conditions S, and V, . 

h h 

An equivalent statement of the least squares formulation involves the 
bilinear form 

(2.9) Bp((\J»,v),(£,w)) 

and the functional 

(2.10) F o(5>w) = f (div w + q£) 

15 fl 

defined for (\|>,v) and (5>w) defined in the finite dimensional subspaces 
S h X ^ of H (fi) X ft ,P (fl). In particular, (^u^) € s h x is the 

minimizer if and only if 


= / (grad - ^)*(grad £ - w) 

a 


+ / r^(div _v + q^)(div w + q£) 

a 


( 2 . 11 ) 


= a11 (5, - } € s h x ^h 
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after selecting a basis for S h x ^h* reduces to a set of algebraic 
equations [9]. As noted in the first section, the chief virtue of this system 
is that it is Hermitian definite. 


3. Approximation and Regularity 

In this section we develop the approximation and regularity results that 
will be needed for the error analysis in the next section. It will be assumed 
throughout this paper that the function q is bounded away from eigenvalues 
of -A (with Dirichlet boundary conditions). Thus (1.1 )— ( 1 .2) has a unique 
solution. Of fundamental importance is the following result due to Kondratiev 
[ 10 ]. 


Theorem 3.1: Let $ satisfy (1.1 )— ( 1.2) on the region ft shown in 

Figure 1-1. Then for t>0, t+l>3>t- + 1 


(3.1) 
and 

(3.2) 


'«,«< c,£, t.e 


< c,£, t.r 


We now turn to approximation. The starting point is a special 
triangulation first introduced by Babuska, Kellogg, and Pitkaranta [7], Three 
conditions are needed in order to obtain the appropriate approximation 
results. The first is the standard angle condition on the individual 
triangles [9]. To describe the second let h be the maximum mesh spacing 
associated with the triangulation. Then, given h and weight factor 3, each 
triangle T must satisfy 
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( 3 * 3 > d(T) = 0(hr g ) 

where 

d(T) = max |x-y|. 
x,y£T 

The third and most crucial rule governs the amount of refinement. 
Given h and 3, the triangles which have the corner as a vertex must satisfy 

1 

(3«4) d(T) £ Ch 1 " 5 . 

For a mesh refined according to these conditions, Pitkaranta [11] proves the 
following. 

Theorem 3.2: Let 0 3 < 1* Let be the space of continuous 

piecewise linear polynomials. Let be defined on Q such that 

(3.5) / r 2 ^ |D 2 i|j| 2 < oo. 

SI 

Then there exists a constant C depending on 3 such that 

(3.6) min {/ |D A («(> - i|> )| 2 + h -1 / r -S | - 4 . | 2 } < ch 2 / r 2P Id 2 41 I 2 . 

n r a ' 

r h h 

Proof : See [11]. 


Equation (3.6) implies that 



- 8 - 


(3.7) 


min 

*h £S h 


'♦ - t.I, < Chlfl, 

hi — £. y a 


We will need a similar result for approximation in other norms. 


Lemma 3.1: For 0 < 8 < 1 


(3.8) min {/ 
\ €S h “ 


* h | 2 + h 2 / r 23 ID 1 ^ 

n 


VI 1 < Ch / 

£1 


: 4S |D 2 *| 2 


Proof: See [12], 


Lemma 3.1 contains two approximation results, namely 


(3.9) 
and 

(3.10) 


>*-Vo< Ch '♦ l 2,2B 


'♦ - *h , l >S i Ch " l ' , 2 , 26 - 


Our analysis will require estimates in dual norms. In particular, 
following will be important: 


(3.11) 






new 


28 , 
r n^ 

1,28 n lnl l,28 


sup / 


The following is an inequality that will be useful in the sequel. 




the 


Lemma 3.2: Ddiv vB 



Proof : Follows from (3.11), Schwarz's inequality, and the inequality 


(3.12) IIV(r 23 n) II 2 < Cllpll 2 2p . 

The final result that will be needed is the analog of the Grid 
Decomposition Property introduced in [6]. This is a condition on ^ , and not 
all finite elements spaces have this property as we will indicate in the 
sequel. The version we will need can be stated as follows. There is a 

number C, 0 < C < 00 and independent of such that for each v^ in 

there are w^z^ in for which 


v, = w, + z, 
h h h 


with 


and 

(3.13) 


div z, = 0 I z, *w, = 0 

h a h h 


lw. I. < Cldiv v, r 
h 0 — h * , g 


4. Error Estimates 

The analysis given here has a structure similar to that found in the 
analysis of mixed methods (see, for example, [6]). The first step is use the 
basic orthogonality property derived from (2.11) to get an estimate in a 
nonstandard norm; in this case it is 
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where B p (»,«) is the bilinear form defined by (2.9). One then uses this to 
derive estimates in the negative norm (3.11), from which one can obtain L 2 
estimates for <p — (J>^. These plus the grid decomposition property yield L 2 
estimates for u - u^. In the sequel we let 

(4.2) 
and 

(4.3) 

In addition throughout we shall assume 

(4.4) 1 - -rr/20 o < B < 1. 


£ = <t> ~ <f> v 


e = 


u - u, . 
— — n 


Lemma 4.1: 


(4.5) | | |(e,e)| | | < Ch{«fll 0}P + II f II x >23 1- 

Proof : It follows from (2.11) that 

(4.6) B p ((e,.e),(iJ; h ,v h )) = 0 
for all (i^v^) € S x $ . Thus 

(4.7) I I I (e,_e) | | | j< inf | | | (<j> - \|> h> u - v^) | | | 

where the inf is taken over all (^jV^) € x V^. Using the approximation 


properties in Section 3 we obtain 
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(4.8) 


e,e)| 1 1 < Ch{nuii 1>g + hii 2>p + nun 2>2e }. 


The inequality (4.5) follows from (4.8) using Theorem 3.1. We now establish 

se' 


an error estimate for div e + qe in the dual norm II * II ^ Q . 


Lemma 4.2: 


(4.9) II div e + qell* < c| | | (e,e) | | | . 

»P 

Proof: For n c W^’ 2 ^ solve 


(4.10) A? + q5 = q in 0 


(4.11) 5=0 on dQ 
with 

(4.12) p = grad 5* 


Observe that by orthogonality and (4.12) 

(4*13) B((e,e),(£ - 5 h ,p - p^) = / r 2 ^ q(div e + qe). 

Therefore, 

(4.14) /r 2B n(div e + qe) _< | | | (e,e) | | | I I I (5 “ 5 h »P - P h )| | | 


< lll(e,e)||| {n? - 5 h » 1 + lip - P h « 1>0 } 
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We choose so that 


(4.15) 
and 

(4.16) 


■ Vl < Chll?ll 2>B < Chlip 3}2p 


1 P - P h 1 1 , (3 — Ch 11 P 11 2 , 2 3 * 


Using regularity (3.1) and (3.2), the estimate (4.9) is obtained by taking 

the sup over n in (4.14) with II n II , „ < 1. 

1 

We now state and prove our two main results. 


Theorem 4.1: Assume that q is bounded on ft. Then 


(4.17) 


I! e H 0 < Ch | | | ( e, e_)J | |. 


Proof: Solve 


(4.18) An + qn = e in ft 


(4.19) n = 0 on 8ft 


for n. For any n^ € we have 

(4.20) Bp((e,£),(0,n - n h )) = Bg((e,£),(0,n)). 

But 

(4.21) 

t 28 

B B ((e,e),(0,n)) = J [grad e»grad n - e*grad r\ + r qn(div e + qe)]. 
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On the other hand 

2 

(4*22) / e = / € (An + qn) = /[-grad e»grad n + qqe]. 

$2 n 

Since 

(4.23) / e*grad q = -/ div eq 

n n 

putting (4.21) into (4.22) gives 

(4.24) J £ = -B ((e,e) , (0,q - q )) + / (1 + r^ q)*(div e + qe). 

fl p 11 SI ~ 

To estimate the right hand side of (4.24) we note that tv can be chosen such 

n 

that 

(4.25) B p ((e,e),(0,q - q h ) < | | |(e,e)| | | | | | (0,q - q h ) | | | 

< ch,, 'T' 2> 3 I I |(e,e)| ||. 

Also 

(4.26) llnll < C II e n _ . < C II £ II n . 

/ ) p U , P U 

Finally, 

(4.27) / (1 + q)n(div e + qe) < clldiv e + qell, Hr 2 ® nil. Oo . 

>p i ,4p 

Moreover, 

(4.28) ,r- 26 "' 1>26 <C< n l 1>0 < 0« e ._ lj0 <C.e. 0 . 


Combining these we obtain (4.17) 
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Corollary: Assume q is bounded on Q. Then 


(4.29) 


11 div £ 11 * A < ch l I |(e>e)| | | . 

> P 


Theorem 4.2: Suppose q is bounded on ft and suppose x 

satisfies the grid decomposition property (3.13). Then 


(4.30) 


llejl 0 1 Ch | | j ( e , e.) | | | + C ^inf llu - ^Hq. 

iih eV h 


Proof : Let 6 S iven * Use (3.13) to decompose as 


follows : 

(4.31) 
where 

(4.32) 
and 

(4.33) 


Uh - Uh = Wh + Zh» 


div = 0, / w^-z^ = 0 


— h h 


IIw^IIq < Clldiv(u h - 11^)11*^. 


Note that for any e x we ^ ave 


(4.34) 0 = Bg((e,£),('|' h ,v h )) = / (grad e - eHgrad <J> h - jv^) 


+ / r^(div £ + qe)(div + q^). 

n a 


Letting ^ = 0 and we have 


^•c = 0. 


(4.35) 
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Thus 

(4.36) l ^I*5h = /Sh<«h ~ 4> = /5h*<“ - 4)’ 

i.e. , 

(4.37) IIz^IIq — “ •Hh’^O* 

But from (4.33) and Lemma 3.2 

(4.38) Iwjjl 0 < Clldiv(u h - 4>"*,3 < Clldiv ell* >g + Clldiv(u - ^)ll* g 

<. Clldiv ejl* g + Cllu - II Q - 

It follows that 


(4.39) 


II ell 


0 i ~ Hh l! 0 + 


*Hh ~ -^h"o < " u “ Vo + "^h“0 + "-Eh “o* 


Combining (4.37)— (4.39) with the above Corollary we obtain (4.30). 


5 . Numerical Results 

In this section we report results of computations which demonstrate the 
weighted least squares method and confirm the analytical results of the 
previous section. All numerical experiments were performed on a VAX 11-750 
computer. Special attention is given to the roles played by mesh refinement 
and the weight. 

All of the examples deal with the Laplace equation 
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(5.1) A* = f. 

We actually solve the equivalent first-order system 


(5.2) 

div u = f 

(5.3) 

u - grad <{> = 


The insensitivity of least squares to type of boundary condition (Dirichlet, 
Neumann, or Mixed) has already been demonstrated [2]. Thus it is sufficient 
for the examples reported here to use the Dirichlet condition 

(5.4) «f> = g on T. 

Consider the L-shaped membrane shown in Figure 5-1. 



( l I, 
2 2 


r 


Figure 5-1: L-shaped Membrane. 
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The re-entrant corner has measure 0 Q = — .■ Thus from [13] we know that 
the solution of (5*1) with homogeneous boundary condition has a singularity 
with leading term 

(5.5) i|> = r 2 ^ 3 sin[-j (0 - tt) ] 

where (r,0) are standard polar coordinates. Therefore, we use rp in (5.5) 
as our test solution. Analysis in Section 4 indicates that optimal rates will 
be assured in the weighted least squares solution by the proper choice of 
weight exponent 0 and correct mesh refinement. The approximating space for 
both <j> and u_ is the space of continuous piecewise linear polynomials. For 
Bq = (3.4) and (4.4) tell us we need 

h . = 0(h 3 ) and ~ < B < 1. 

min J 

Symmetry allows us to solve on the region shown in Figure 5-2, with a 
tangency condition imposed in the line of symmetry, as shown. 



Figure 5—2: Computational Region for L-shaped Membrane. 
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Figures 5-3 and 5-4 display two of the finite element grids used in the 
numerical experiments. Note that each triangulation is constructed by 
subdividing the basic criss-cross grid so that every element is an isosceles 
right triangle. This type of refinement was chosen instead of the coordinate 
stretching from a uniform mesh since the latter contained some elements with 
large aspect ratios. To assure accurate determination of convergence rates, 
the meshes were constructed so that the number of points N varied with the 
maximum mesh spacing h according to the relation 

N = 0(h“ 2 ). 



Figure 5-3: Refined Mesh, h =1/2 
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Figure 5-4: Refined Mesh, h = V 4 . 

The L 2 errors in u and <{> are displayed in Figure 5-5. Results from 
the weighted least squares scheme on a refined mesh are contrasted with those 
using the standard least squares scheme on a uniform grid. 

We also applied the weighted least squares scheme to Laplace's equation 
on a square region with a crack, illustrated in Figure 5-6. This model is 
characteristic of such physical problems as torsion of a cracked beam and flow 
over a very thin airfoil. 

It must be noted that the analytical results do not hold for this problem 
because the cone condition used in the regularity result is not satisfied. 
However, the results hold if the crack is replaced by a re-entrant corner with 
measure 9 q = 2ir - e. For this problem, (3.4) and (4.4) indicate that the 
mesh refinement and weight parameters must satisfy 

h min = 0(h4) and f < & < !• 
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Figure 5-5: 1.2 error in weighted least squares approximation to $ (a) 

and u_ = 7<(i (b) contrasted to standard (unweighted) least 

squares approximation to <(> (c) and u (d) . 
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Figure 5—6: Cracked Square. 

For the test solution we again use the leading term in the singularity 
which is 

(5.6) = r^^ sin[-j (® “ 7T )]* 

As before, by symmetry we will solve only on the region shown in Figure 

5-7. 



Figure 5-7: Computational Region for Crack Problem. 
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The L£ errors in the weighted least squares solution to the crack 
problem are displayed in Figure 5-8. 

The following conclusions can be drawn from these computations. The 
accuracy lost in the least squares solution due to a corner-type singularity 
is fully restored by the use of a weighted scheme on a refined mesh. Specific 
criteria, which depend on the measure of the corner, have been developed to 
determine the correct weight and order of refinement. 

Moreover, the weighted scheme inherits all the advantages associated with 
least squares. Second-order accuracy is achieved using the same finite 
element spaces for <J> and jj. The associated matrix system is always 
symmetric and positive definite, allowing solutions by standard iterative 
techniques. The essential boundary conditions can be included in the 
variational principle instead of being imposed directly on the approximating 
finite dimensional space. 
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